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INTRODUCTION 


The family of beta distribution with density 


/ 


f(x|a,b,a,B)~ r(a)T(B) 


r(a+B) 



with range parameters a<x<b and shape parameters a,B>0 is an extremely rich 
family of distributions. It includes J shaped, U shaped and unimodel distri- 
butions within finite bounds. In spite of such versatility and applications 
in physical and life sciences, tables and computer programs for probabilities 
for the beta distribution are generally not available. One might suspect 
that this is because the standard beta distribution (a=0, b=l ) probabilities 
and fractile points can be computed from F distribution tables which are 
generally available. However, this is not practical for even a moderate num- 
ber of probabilities or fractile points, and is limited to integral values 
of a and B. In practice beta probabilities and fractile points are often 
needed. The enclosed computer program meets this need. 

The program presented here computes the probability that a beta random 
variable, X with shape parameters a and B» on the range a to b will fall 
within an interval (x x , x 2 ) a£Xj£x 2 £b. Also for any given probability p the 
program will compute the fractile x p =F -1 (p)* Additionally routines to cal- 
culate the density function f(x) and the gamma function F(x) are presented. 
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a/79- Z/7^7 


SYMBOLS 


r (x) 

X, V 
a 

x, y 

F 

f, b 

r 

t 

b, a 

a> 8 
n 15 n 2 

X 2 

l 


Gamma Function 

. > “ -t jc-l 

F(x)= Set dt 
0 

Random Variables 

Real Number 

Real Numbers 

Distribution Function 

Density Functions, or Real Numbers 

Root of a Real Valued Function 

Dummy Variable for Integration 

Range Variables for Beta Distribution 

Shape Parameters for Beta Distribution 

Degrees of Freedom for F-distribution 

Chi Squared Test Value 

Summation Notation 


FORTRAN FUNCTIONS 


BDIST 

BFRAC 

F 

G 


Beta Distribution Function 
Beta Fractile 
Beta Density Function 
Gamma Function 


PROCEDURE 


The computation of probabilities and fractiles of beta distribution is 
based on the following procedure. The general beta density for random var- 
iable x is given by 

f(x| “’ 6)= r(Sr(s) (f|) {&£f (b^i) 

The range of variable x is interval (a,b) and the shape parameters are a and 
3>0. The transformation y=(x-a)/(b-a) yields standard beta density with 
range (0,1) and density 

which takes different shapes (see figure 1) depending, on the values of a and 
3. These include J-shaped, U-shaped and symmetric and skewed unimodel dist- 
ributions. 



Figure 1. p(y) = y ° 1 ^ ~ 1 ’ the probability 

density for various values of a and p . 
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The beta probabilities are found by numerically integrating the beta 
density function. The technique for this is Gauss-Legendre quadrature 
(Ref. 1). The integral is translated to the interval (-1,1) and the Legendre 
polynomial of order 15 is used. 

Having determined the integral of the beta density function, the p-th 
fractile point Xp of the beta distribution is obtained by solving the inte- 
gral equation 


F(x )= / P f(t)dt = p 

-00 


The distribution function F(x) is always a monotone-increasing function. 
As such, virtually any iterative algorithm will converge to the solution of 
the equation F(x)-p=0. The technique used in the programs presented here is 
the method of bisection (interval halving) (Ref. 2). The procedure termin- 
ates at x whenever | F ( Xp ) -p | <1 0~ 6 . 

The gamma function F(x) for argument x is evaluated by using sterling's 
approximation (Ref. 3) 

(*) T(x) = (2ir/x) exp (xlnx-xg(x)) 

where 

g(x) = l-l/12x 2 + l/360x 4 - l/1260x 6 + 1/1 680x e 

us i ng Wy+mI / \ 

(**) r ( x > =iJ ^0 x { } = x(x+l)--*(x+n-l) 

n can be selected so that x+n>10 then T(x+n) computed by the series expan- 
sion, (*) and T(x) found by (**). From the calculation of T(x) the beta den- 
sity function is easily determined. Appendix A lists the FORTRAN routines. 

OUTPUTS 

The four outputs available are: 

P(x!<X<x 2 ) where X is a beta random variable 

x =F 1 (p), where F is a distribution function of a beta random 

r 

variable 

f(x|a,b,a,3), the density function of a beta random variable 
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T(x), the value of the gamma function for positive argument, x. 
These outputs and the parameters for each FORTRAN function are described 


below. 


Beta Distribution Function 

The distribution function value F(x) for the beta distribution on the 
range (a,b,) with shape parameters a and $ is given by 

BPR0B(0. , x, ALPHA, BETA, A, B). 

The probability a beta random variable with parameters a,b,a, and B takes a 
value between Xj and x 2 is given by: 

BPR0B(Xj , x 2 , ALPHA, BETA, A, B). 


Example : The probability that a beta random variable with range (1., 3.5) 

and shape parameters a=4.1 3=16.2 takes a value less than 1.75 is 

BPR0B(1 . , 1.75, 4.1, 16.2, 1., 3.5) = .864 . 

Beta Distribution Fractile Points 

For a given p, the lOOp-th fractile point from the beta distribution 
with parameters a,b,a, and 6 is given by: 

BFRAC(P , ALPHA, BETA, A, B) 

Example 1 : The 75-th percentile of the beta distribution on the range 240 to 

1400, with shape parameters ALPHA=14.2, BETA=34.7 is 

BFRAC( . 75, 14.2, 34.7, 240., 1400.)= 559.307 

Example 2 : Appendix B presents a FORTRAN program which uses BFRAC to gener- 

ate ten equi probability points of a beta distribution on the range (0,, 1.) 
with shape parameters ALPHA=4., BETA=16. Thirty beta ramdom variables from 
this distribution are generated and a x 2- test is performed, at the 5 per cent 
level, to test the null hypothesis that these observations are, indeed, from 
the standard (i.e. range(0., 1.)) beta distribution with shape parameters 
ALPHA=4., BETA=16. Not surprisingly, the x 2- test statistic accepts the null 
hypothesis. 


Beta Density Function 

The probability density function f(x|a,3»a,b) for a beta random variable 
on the range (a,b) with shape parameters a and B is given by: 

F ( X , C, ALPHA, BETA), 


where C= T(ALPHA+BETA)/(r(ALPHA) T(BETA)) must be passed to the function. 
C may be computed using the gamma function 
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C= GAMMA (ALPHA+BETA)/ ( GAMMA (ALPHA) *GAMMA (BETA)) 

Example : The beta density function on the range (.67 to 12.2) with shape 

parameters ALPHA=1.2, BETA-3.5 at x=5.2 is given by: 

C=GAMMA(4.7)/(GAMMA(1.2)*GAMMA(3.5)), F(5.2, C, 1.2, 3.5)= 0.1045. 
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APPENDIX A 


This appendix presents the computer proarams used for computinq the 
beta-density probabilities and fractile points. 
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FUNCTION FFRAC(X, ALPHA, 9FTA,A,B) 

Q+ ¥ ¥ ¥ 

C**** FOUTINE TC DETERMINE FRACTILE POINTS OF 

C **** THE BETA DISTRIBUTION* WITH PARAMETERS 

C**** ALPHA ANO BETA 

Q ¥* ¥ ¥ 

C**** I N VERSE DISTRIBUTION VALUE FOUND EY 

C**** HAL VI NG- THE -INTER VAL TECHNIQUE* . 

Q**** 

C •*** BY BROWNLCW, SOC/ISI. 

C*** * 

c**** MAKE 5 URE PERCENTAGE POINT IS BETWEEN 0 AND 

IFC X ,Ll. 1. • ANO. X . GE. 0. ) GO TO 2 

PRINT l f X 

1 FORMAT (•• X = Ell. 3 , '* FOR BETA OIST*N, SO FRACTILE POINT COULD N 
-OT BE FOUND*') 

£ ¥ ¥ ¥ ¥ 

bffa: = o. 

RETURN 

2 KOUNT = 0 

TOF = 1. - 1.E-5L 
30T = 1 .E-5D 

£¥ ¥ ¥ ¥ 

3 aFSA: = (TOF ♦ BOT ) / 2 • 

Pi = BFR06U BFRAC, ALPHA, BETAtC* ,1. ) - X 

£ ¥ ¥ ¥ ¥ 

C**** IF THC- INTEGRAL OP THE BETA FUNCTION FROM ZERO TO 

C**** BFF AC DIFFERS FROM X BY LESS THAN l.E-9 

C**** WE HAVE THE ERACTILE POINT.. 

C*** * 

IF ( A B S (FI) • L T • l.E-9) GO TO 4 

C ¥ ¥ ¥ ♦ 

C *♦** OTHERWISE RESET THE UPPER OR LOWER INTERVAL 

£**•* VALUE AND REPEAT THE PROCEOURE.. 

[»¥¥¥ 

IF C r l .GT . *. ) TOP = BFRAC 
IF ( -1 .LT • 5. ) BOT = BFRAC 

Q¥ ¥ ¥ ¥ 

C**** IN ANY CASE NEVER DO MORE THAN 20 ITERATIONS 

C**** TO FIND A VALUE ON (0.,1.) TO WITHIN 2**(-19). 

Q ¥ ¥ ¥ ¥ 

KOUNT = KOUNT ♦ i 
I F ( KOUNT .GT. 20 GO TO 4 

£¥ ¥ ¥ ¥ 

GO T3 3 

Q¥ ¥ ¥ ¥ 

C**** REDEFINE FRACTILE POINT FCR DISTRIBUTION ON THE 

C **** RANGE A TO B.. 

£¥¥ ¥ ¥ 

4 BFPAC = A ♦ BFRAC* (B-A ) 

C¥¥¥¥ 

C ¥ ¥ ¥ ¥ 

RETURN 

ENO 
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C*** * 

c**** 

C** 4 ** 
C**** 

c**** 

C**** 

c**** 

(;♦♦♦♦ 

C**** 

Q4HMMF 

f * 4 * ¥ 

C**** 

C**** 

Q* ¥ ¥ ¥ 
£ ¥ ¥ ¥ ¥ 


1 


c **** 

? 

3 


C ¥ ♦ ¥ ¥ 

4 

5 


C** * * 

6 

7 


C**** 

8 

C**** 
C**** 
C**** 
C*** * 


FUNCTION EFROB (XI, X2, ALPHA , BET A, A. 9) 

FOUTINE TO INTEGRATE BETA DENSITY FUNCTION 
DEFINED ON ( A 1 3 ) WITH PARAMETERS ALPHA* AND 
BETA, LIMITS OF INTEGRATION* Xi TO X2. 

BY 3ROWNLCW, SQC/ISI, 1/79 


EXTERNAL F 


F IS THE BETA DENSITY FUNCTION, 

F (XtCt ALPHA f BETA) = CM ( (X-A) / (0-A >)*M ALPHA-1 .) * 

( (B-X)/ (B-A)** ( BETA-1. ) )/ (B-A) 


C = GAMMA (ALPHAS BET A) / ( GAMMA (ALPHA) ♦GAMMA (BETA ) ) 


UNDER THE TRANSFORMATION Y = (X-A)/(B-A), 

Y HmS STANDARD DENSITY, 

F ( Y ) = CM ( Y * * (ALPHA -1 •) * ( 1 • - Y ) * * (BETA-1 ) ) 


IFf A . LT, B ) GO TO 2 
PRINT i,A,6 

FORMAT (* IN ESTIMATING BETA FRO 3A 31 L IT Y , THE LOWER BO UND , * , El 2 . 4 , 
/* rs NOT GREATER THAN THE UPPER 30 LNC* , * , El 2. 4) 

BP p 03 = 'i, 

RET UR U 

IF (XI .LT, X2 ) GO TO 4 
PRINT 3 , XI , X2 

FORMAT <♦ IN ESTIMATING BE T A PROBABILITY, THE LOWER L I MI T , * , El 2 . 4 , 
/* IS NOT GREATER THAN THE UPPER LIMIT ,* ,E12, 4) 

3PP0B - ). 

RETURN 


IF ( <1 . GL . A) GO 1C 6 
PRINT f ,Xj , A 

FOFMA T ( * IN ESTIMATING THE BETA PROBABILITY, THE LOWER LIMIT,*, 
E12.4,/,* IS LESS THAN THE LOWER BOUND OF THE DISTRIBUTION,* 

El ?• 4 , * 7 HE LOWER BOUND WILL BE USEO*) 

XI = A 


IF(X2 . LE, 8) GO TO 8 
PRINT 7 , X2 , B 

FOP MA T ( * IN F STIMATING THE 3FT A PROBABILITY, THE UPPER LIMITS 
E12.4,/,* IS GREATER THAN THE UPPER BOUND OF THE DISTRIBUTION* 
El 2 • 4 , * THE UPPER BOUNC WILL BE USED,*) 

X2 = 3 

Y 1 = ( X 1 -A ) / (8-A) 

Y2 = ( X 2- A ) / ( B-A ) 

GAUSSLQ DOES GAUSS I A N-L EGENDRE QU AORITURE TO 
ESTIMATE THE INTEGRAL OF THE BETA OENSITY FUNCTION. 

CALL G^USSLO(P, ALPHA, BETA, Yl,Y2t ANS ) 

BPF03 = AKS 


RETURN 

END 


SUBROUTINE GAUSSLQ (F t ALF HA , 3ETA , XI f X 2 , ANS) 

C**** 

C * ♦♦♦ ROUTINE TC HO GAUSS I A N- LEGENDRE OUAORITURE 

C ♦♦♦♦ TO INTEGRATE BETA DENSITY FUNCTION FROM XI 

TO X2 

C**** F IS THE STANOARO BETA DENSITY FUNCTION.*. 

C**** F (X»C» ALF HA * BET A ) WHERE C IS NORMALIZATION 

C ♦♦♦♦ CONSTANT, C= GAMMA (AL PHA ♦■BETA )/ (GAMMA (ALP HA) 

♦ GA MM A ( BE TA ) ) 

C ♦*♦♦ FY 3ROWNLCW, SOC/ISI, 3/79 

Q# ¥ * * 

01 MF MS I ON FOOT S ( 3) ,W(3) 

DOUBLE PRECISION ROOTS , W , ANS WE R 

C = GAMMA (AL p HA f BET A )/( GAMMA ( ALF HA ) ♦GAMMA ( 8E TA ) ) 

DATA RCOTl/ 


- C.JOl 


, C.2C119 

4 J939 

97435 OC , • 

39415 

1 3 47H 

7756300 , 

- C.57C37 

21726 

J 36390.*, .72441 

77313 

6C1700G , • 

34320 

65 334 

104270b , 

- (.93727 

3 3 92 4 

7ooOi , . 98799 

2 5131 

204 35 01 / 




C* * # * 








DATA W/ 

- C.2P257 

3241 9 

255610C, . 19343 

1 4353 

2711100,. 

13616 

10 001 

1556200 , 

-■ C.i6i26 

9 2j{ 3 

lb9940o, .13997 

C 6779 

2 61 540T , • 

10 715 

92 2u 4 

67172DL , 

- i.37f3<b 

6)474 

3 31 v 3D J, . 0 30 75 

32419 

96117DC / 





0 ¥ ¥ ¥ ¥ 

0 ¥ ¥ ¥ ¥ 

C"" TRANSLATE THE INTEF/AL (Xl,X2) TO (-1,1) 

0 ¥ ¥ ¥ ¥ 

ANSWER = mi) * r ( (XifX2) /2. ,C, ALPHA, BETA) 

0 ¥ ¥ ¥ ¥ 

C ¥ ¥ ¥ ¥ 

JO 1 1=2,8 

0 ¥ ¥ ¥ ¥ 

ANSWEP = ANS WEP ♦ W ( I ) ♦ F ( ( ( X2-X1 ) ♦ ROOTS ( I ) ♦ ( XI «• X2 ) ) / 2. ,C , 

ALF HA , BET A ) 

0 ¥ ¥ ¥ ¥ 

ANSWER^ AN! we° «* WII)» F ( ( ( X2-X1 ) ♦ ( - 1. ) ♦ROOTS ( I ) ♦ < XI + X2) ) /2. , C, 

ALF HA , BET A ) 

1 CONTINUE 

0 ¥¥ ¥ ¥ 

0 ¥ ¥ ¥ ¥ 

ANS= A NS W = R ♦ ( X 2 -X 1 ) /2. 

0 ¥ ¥ ¥ ¥ 

RETURN 

ENT 
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FUNCTION F<X»C»ALPHA,8ETA) 


C**»* 

C**** F IS THE BETA DENSITY FUNCTION... 

C 

C**"* C FASSEO IN MUST HAVE THE VALUE 

C**** C= GAMMA(ALPHA-*-3ETA)/(GAMMA (ALP HA) 4 GAMMA ( 3ET A) ) 

0 * ¥ 

C**** 8 Y 8R0WNLCW, SOC/ISI, 1/18/79 

0 

F = 3* ( X* * (ALFHA-1.)* ( 1* -X) ** ( 9ETA -1 • ) ) 

(;♦♦♦♦ 

RETURN 

ENG 
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FUNCTION GANNA (X) 


C+*** 

C 

C**'* 

c**»* 

C **** 

C***» 

C*** * 

£¥*♦♦ 

C**»* 

C*** * 

C***» 

C*MMM> 

£«##« 

Q* + + + 

c*+* + 

Q* + + * 

A = X 
P = 1. 

C**** 

1 IF ( A .GT. 10, ) GO TO 2 
P = P»A 

A = A ♦ 1. 

GO TO 1 

C** + * 

£*¥*♦ 

2 CONTINUE 

C**** 

E- 1-<1./<12.*A*A) )*i./( 36Q.*A*A*A*A> - l./<1260.* 
. A*A* A* A* A*A) 4- l./(1650.*A*A*A*A*A*A*A*A) 

(;♦♦** 

OATA PI/3. 1415926538/ 

GAMMA = ( SORT (2 • * F I / A ) * EXP< A* < AL OG ( A ) -E ) ) ) /P 

C**** 

c**** 

RETURN 

END 


WRITTEN BY BROWNLOW.SDC/ISI . JUNE 1578 

GAMMA FUNCTION. REAL ARGUMENT MUST BE 
GREATER THAN ZERO. 

G ( X ) = INTEGRAL ( EXP ( -T ) *T** ( X-l) DT 
E tfALUATEO BY EXPANSION OF STIRl ING f S 
F CRMULA 


GAMMA (X ) = GAMMA(XfN) /(X* (X*l) .♦ (X*N) ) 
WHERE X 4- N IS CHOSEN TO GET MAXIMUM 
ACCURACY FROM THE ALGORITHM. SEE 
HENRICI. “COMPUTATIONAL ANALYSIS 
WITH THE HP-25 POCKET CALCULATOR" 
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APPENDIX B 


This appendix presents the FORTRAN program which generates the ten equi- 
probability points of the standard beta distribution with shape parameters 
a=4, 8=16. Additionally thirty random variables from this distribution are 
generated and a x 2 goodness-of-fit test is performed. 


PROGRAM MA IN (C UT PUT) 




ROUTINE. TC FIND ORDINATE VALUES' 

c + *** 


OF THE BETA DISTRIBUTION WHICH 

Q + + ** 


DIVIDE THE (0., i.) INTFRVAL 

Q * ¥ * ¥ 


INTO TEN SHUNTER VALS OF FOIAL 



FRO 3A8IL ITy . . 

c **** 


BETA DISTRIBUTION WITH ALPHA = h. 



BEI A = 16. 

£¥¥¥¥ 

c**** 


6 V 3RGWNLCW, SOC/ISI. 

c ++* + 




DO 2 1=1 

» q 


X = FLOAT C I ) / 1 0 • 


PRINT i. 

BFRAL (X,4. ,16. .0. .1.) 

i 

FORMAT (* 

*F10. 5) 

2 

continue 


c **** 


AND NOW GENERATE 3G SA-IPLIS 

c **** 


FROM A STANDARD BETA DISTRIBUTION 

c **** 


WITH FARAFTCRS ALPHA = 4., 

£ ¥¥ ¥ ¥ 


AND BETA - 16. 


no 4 1=1 

,30 


PRINT 2, 

P E TAGS N (4,16,0. ,1.) 

3 

FORMAT ( * 

* F5.3 ) 

4 

CONTINUE 



END 
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FUNCTION 8ETAGFN (I ALPHA, I 8FT A* A » 9) 

£¥ ¥ ¥ ¥ 

c*** * 

£¥ ¥ ¥ ¥ 

c»*»* 

c»»»* 

£¥ ¥ ¥ ¥ 

C ¥ ¥ ¥ ¥ 

£¥ ¥ ¥ ¥ 

Q¥¥ ¥ ¥ 

C ¥ ¥ ¥ ¥ 

£¥¥ ¥ ¥ 

G1 

DO 
G1 

1 CONTINUE 

C¥¥¥ ¥ 

G2 = C. 

DO 2 1 = 1 ♦ 1 9fc' T A 

G2 = G 2 - A L 0 0- ( RANF(U) ) 

2 CONTINUE 

C ¥ ¥ ¥ ¥ 

BEDAGEN = A * (G1/(G1 * G2))¥(D-A) 

C¥¥¥¥ 

C ¥ ¥ ¥ ¥ NOTE THAT Gi AND G c ARE* IN 

C ¥ ¥ ¥ ¥ FACT, GAMMA FANDOM VARIABLES. 

Q¥¥ ¥ ¥ 

St TURN 
£NC 


FUNCTION TO GENERATE A t'ANDCM 
V'ARAIBLE FROM A 8-TA 
DISTRIBUTION CN T ML RANGE 
(A, 9) WITH SHAPE PARAMETERS 

IALPHA, AND I q ETA ♦ 

METHOD IS FROM HASHNG3 AND 
FFmCOCK, "STATISTICAL DIST- 
RIBUTIONS". PAGE 72. 

= 0 . 

1 1=1, IALPHA 
= Gl - A L OG ( RANF(O) ) 


The equiprobabil ity points for this distribution which partition the range 
(0,1) into ten intervals are: .09514, .12334, .14652, .16817, .18989, 

.21297, .23907, .27131, .31895, 1.00000. 

The thirty beta random variables generated are: .096, .238, .147, .230, 

.234, .136, .202, .074, .338, .153, .216, .167, .370, .350, .167, .192, .171, 
.229, .144, .056, .257, .201, .137, .235, .181, .339, .156, .122, .235, .137. 


15 


Given the equi -probability points and the random sample found by using 
the computer programs given in this paper, the null hypothesis that the dis- 
tribution from which these random variables are drawn is beta on (0,1) with 
shape parameters a=4, 8=16. may be tested. Specifically, 


Class 


Observed Freq. 


Expected 


0.00000-0.09514 2 
0.09514-0.12334 2 
0.12334-0.14652 4 
0.14652-0.16817 5 
0.16817-0.18989 2 
0.18989-0.21297 3 
0.21297-0.23907 6 
0.23907-0.27131 1 
0.27131-0.31859 0 
0.31859-1.00000 4 


3 

3 

3 

3 

3 

3 

3 

3 

3 

3 


Each interval has equal probability (in this case 0.10) so with a random sam- 
ple of size thirty we expect three observations in each interval. 


The test statistic is given by 


x 2 = 


2° (0,-E.) 2 /E. 
i=l 1 1 1 


t h 

where 0- is the observed number of occurances in the i cell; E- the expect- 
ed number. 

1 


hence X 2= J (1 2+ l 2+ l 2 +2 2 +* • *+l 2 ) 
= T (31 )=1 0. 333 


and x 2 ’ (9)=16.919 

.05 


thus we accept H_: the distribution from which the sample was drawn is 

6 (4.16) 
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